#How did the prescription of the top 10 most prescribed antibiotics in 2019 change due to lockdown
#How did covid-19 pandemic and lockdown affect the prescription trends of the most common antibiotics in Scotland. (Would it be better to focus on the main lockdown (May 2020) or also compare between lockdowns?)
First I will load the appropriate libraries. Then, load, clean and filter the data.
First, I loaded the datasets and created a function to concisely use the same functions on different datasets. I learned how to do this from the website https://www.dataquest.io/blog/write-functions-in-r/. In this function I filtered for antibiotics by only looking at prescriptions with the BNFItemCode starting with 0501. Then I filtered out NAs, cleaned up the names and calculated the number of prescriptions. I then merged the datasets
may2020 <- read_csv(here("data","pitc202005.csv"))
may2024 <- read_csv(here("data","pitc202405.csv"))
may2019 <- read_csv(here("data","pitc201905.csv"))
may2019 <- may2019 %>%
rename("HBT" = HBT2014)
process_antibiotics <- function(data) {
data %>%
filter(str_detect(BNFItemCode, "^0501"), !is.na(BNFItemDescription)) %>%
mutate(BNFItemDescription = str_replace(BNFItemDescription, "_", " "),
Antibiotic = str_extract(BNFItemDescription, "^[^ ]+")) #Extracts the first word, which is the antibiotic name
}
#I filtered for prescriptions with an item code beginning with 0501 because prescriptions with that item code are antibiotics (Open Prescribing, 2019).
may2019_processed <- process_antibiotics(may2019)
may2020_processed <- process_antibiotics(may2020)
may2024_processed <- process_antibiotics(may2024)
#I split this into 2 functions as later I want to group by HBT
antibiotics_by_name <- function(data2) {
data2 %>%
group_by(Antibiotic) %>%
summarise(quantity_sum = sum(PaidQuantity), .groups = "drop") %>%
arrange(desc(quantity_sum))
}
#I filtered for prescriptions with an item code beginning with 0501 because prescriptions with that item code are antibiotics (Open Prescribing, 2019).
may2019_processed <- antibiotics_by_name(may2019_processed)
may2020_processed <- antibiotics_by_name(may2020_processed)
may2024_processed <- antibiotics_by_name(may2024_processed)
I then prepared the data for joining by renaming Fluclox to Flucloxacillin. I then joined the data, and tidied it
may2020_processed <- may2020_processed %>%
mutate(Antibiotic = ifelse(Antibiotic == "FLUCLOX", "FLUCLOXACILLIN", Antibiotic))
may2019_processed <- may2019_processed %>%
mutate(Antibiotic = ifelse(Antibiotic == "FLUCLOX", "FLUCLOXACILLIN", Antibiotic))
antibiotics1920 <- full_join(may2019_processed, may2020_processed, by = "Antibiotic")
antibiotics192024 <- full_join(antibiotics1920, may2024_processed, by = "Antibiotic")
antibiotics192024 <- antibiotics192024 %>%
rename(
"May 2019" = quantity_sum.x,
"May 2020" = quantity_sum.y,
"May 2024" = quantity_sum) %>%
slice(1:10)
antibiotics192024 <- antibiotics192024 %>%
mutate(Antibiotic = str_to_sentence(Antibiotic)) #I changed all the antibiotic names to only have the first letter capitalised. I learnt this from the R help page
I created a stacked bar graph to visually compare prescriptions
pivoted_antibiotics <- antibiotics192024 %>%
pivot_longer(
cols = c("May 2019", "May 2020", "May 2024"),
names_to = "Year",
values_to = "Prescription_Number")
options(scipen=999)
antibiotics_graph <- pivoted_antibiotics %>%
ggplot(aes(x = Prescription_Number, y = Year, fill = Antibiotic)) +
geom_bar(stat = "identity") +
labs(title = "Prescriptions Sold for Each Antibiotic by Year") +
scale_fill_brewer(palette = "Paired") +
theme_minimal()
antibiotics_graph <- antibiotics_graph %>%
ggplotly(tooltip = c("Prescription_Number", "Antibiotic"))
antibiotics_graph
plotly - interactive graphs in R
This data visulaisation shows that the overall number of prescirptions sold decreased during lockdown and was slightly lower before lockdown than after. It also shows that the top 10 most prescribed antibiotics were not the same in every year.
I then created a gt() table to show the number of prescriptions and calculate the total prescriptions per year
all_antibiotics_table <- antibiotics192024 %>%
summarise(across(starts_with("May"), sum, na.rm = TRUE)) %>%
mutate(Antibiotic = "Total") %>%
bind_rows(antibiotics192024, .)
all_antibiotics_table <- all_antibiotics_table %>%
gt() %>%
cols_label(Antibiotic = "Antibiotic") %>%
cols_align(align = "center",
columns = everything()) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels(columns = everything())
) %>%
tab_header(
title = "Affect of the Covid19 pandemic on the prescription of common antibiotics",
subtitle = "Looking at the change of prescriptions for the top 10 most common antibiotics from before, during and after the longest Covid19 lockdown") %>%
tab_spanner(
label = "Number of Prescriptions",
columns = c("May 2019", "May 2020", "May 2024"))
all_antibiotics_table
| Affect of the Covid19 pandemic on the prescription of common antibiotics | |||
| Looking at the change of prescriptions for the top 10 most common antibiotics from before, during and after the longest Covid19 lockdown | |||
| Antibiotic |
Number of Prescriptions
|
||
|---|---|---|---|
| May 2019 | May 2020 | May 2024 | |
| Amoxicillin | 2875868 | 1778491.0 | 3751932.0 |
| Phenoxymethylpenicillin | 1892108 | 968754.0 | 2616743.0 |
| Flucloxacillin | 1612959 | 1475567.0 | 1733244.0 |
| Trimethoprim | 782422 | 745904.0 | 696670.0 |
| Erythromycin | 661545 | 443204.0 | 582068.0 |
| Lymecycline | 623567 | 493844.0 | 538300.0 |
| Oxytetracycline | 572089 | 428371.0 | 277441.0 |
| Doxycycline | 447389 | 387163.0 | 616170.0 |
| Co-amoxiclav | 382770 | 349958.0 | 466694.0 |
| Clarithromycin | 377144 | 230282.1 | 493464.5 |
| Total | 10227861 | 7301538.1 | 11772726.5 |
I then created a faceted geospatial map to show which areas of Scotland were prescribed more antibiotics than other areas, and how this changed over covid
#Load data
NHS_healthboards <- st_read(here("data", "SG_NHS_HealthBoards_2019.shp"))
## Reading layer `SG_NHS_HealthBoards_2019' from data source
## `C:\Users\franc\OneDrive\Documents\Data Science\Assessment\B207252\data\SG_NHS_HealthBoards_2019.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 5512.998 ymin: 530250.8 xmax: 470332 ymax: 1220302
## Projected CRS: OSGB36 / British National Grid
#Prep for joining
NHS_healthboards <- NHS_healthboards %>%
rename("HBT" = HBCode)
#Summarise antibiotic data by HBT code to give a total per area
antibiotics_by_HBT <- function(data_frame) {
data_frame %>%
group_by(HBT) %>%
summarise(quantity_sum = sum(PaidQuantity), .groups = "drop") %>%
arrange(desc(quantity_sum))
}
#Apply the function
HBTmay2019_processed <- antibiotics_by_HBT(may2019)
HBTmay2020_processed <- antibiotics_by_HBT(may2020)
HBTmay2024_processed <- antibiotics_by_HBT(may2024)
#Join antibiotic data
HBTantibiotics1920 <- full_join(HBTmay2019_processed, HBTmay2020_processed, by = "HBT")
HBTantibiotics192024 <- full_join(HBTantibiotics1920, HBTmay2024_processed, by = "HBT")
#Rename and remove NAs
HBTantibiotics192024 <- HBTantibiotics192024 %>%
rename(
"May 2019" = quantity_sum.x,
"May 2020" = quantity_sum.y,
"May 2024" = quantity_sum) %>%
filter(!if_any(everything(), is.na))
#Pivot to make faceting easier
HBTpivoted_antibiotics <- HBTantibiotics192024 %>%
pivot_longer(
cols = c("May 2019", "May 2020", "May 2024"),
names_to = "Year",
values_to = "Prescription_Number")
#Join with healthboard data
Healthboards_and_antibiotics <- NHS_healthboards %>%
full_join(HBTpivoted_antibiotics, by = "HBT")
#Geom plot
Healthboards_and_antibiotics_plot <- Healthboards_and_antibiotics %>%
ggplot(aes(fill = Prescription_Number)) +
geom_sf(color = "white", linewidth = 0.1) +
scale_fill_viridis_c(name = "Total No. of Prescriptions",
labels = scales::comma) +
labs(title = "Total Number of Antibiotic Prescriptions in Areas of Scotland from 2019 to 2024") +
facet_wrap(~ Year) +
annotation_scale(
location = "tl"
) +
annotation_north_arrow(
location = "tl",
pad_y = unit(0.5, "in"),
style = north_arrow_nautical(
fill = c("grey40", "white"),
line_col = "grey20"
)
)
geog_graph <- Healthboards_and_antibiotics_plot %>%
ggplotly()